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Riccati differential equations have played important roles in the theory and 
practice of control systems engineering. Our goal in this paper is to propose a 
new multistage successive approximation method for solving Riccati 
differential equations. The multistage successive approximation method is 
derived from an existing piecewise variational iteration method for solving 
Riccati differential equations. The multistage successive approximation 
method is simpler in terms of computing implementation in comparison with 
the existing piecewise variational iteration method. Computational tests show 
that the order of accuracy of the multistage successive approximation method 
can be made higher by simply taking more number of successive iterations in 
the multistage evolution. Furthermore, taking small size of each subinterval 
and taking large number of iterations in the multistage evolution lead that our 
proposed method produces small error and becomes high order accurate. 
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1. INTRODUCTION 

Mathematical modelling and simulation have been applied extensively in the areas of 
telecommunication [1], [2], computing [3]-[5], electronics [6], [7] and control [8], [9]. One of important 
models in these areas is the Riccati differential equation. Riccati differential equations occur in control 
systems engineering. Control systems engineering itself has played important roles in electrical engineering 
and related areas [10]. Riccati differential equations are quadratic with respect to the unknown function. 
These equations are named after a Venetian mathematician, Jacopo Francesco Riccati (1676-1754) [11]. The 
term Riccati equation is also used to refer to an analogous matrix equation occuring in quadratic control 
problems. The non-dynamic steady-state version of them is referred to as the algebraic Riccati equation. 

Riccati equations and their properties are applied in recent publications [12]-[15]. A number of 
studies relating to Riccati equations are also reported in the literature [16]-[21]. Due to the importance of 
Riccati equations, we focus on proposing a computing method for solving Riccati differential equations. 

Amongst available methods in the literature, successive approximation and variational iteration 
methods are able to provide accurate solutions near the initial point of the domain. Successive approximation 
methods are successful in solving various problems [22]-[23]. Variational iteration methods are also powerful 
in solving a wide variety of mathematical models [24]-[27]. Interestingly, these two methods (successive 
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approximation and variational iteration methods) are identical for a particular setting, as reported by 
Jafari [28]. 

Our contribution (our goal) in this paper is to propose a new multistage successive approximation 
method for solving Riccati differential equations. We recall an existing piecewise variational iteration 
method due to Geng, et al. [29]. We modify the method of Geng, et al. [29] to obtain the multistage 
successive approximation method. The modification leads to a simpler implementation of the resulting 
method in terms of computation. Riccati differential equations have the following form 


d 
= = p(x) +a )yx) + roy) (1) 
on 0 < x < X, with initial condition 


y(0) = yo (2) 


where x is the free variable, y(x) is the unknown function dependent on x, X is a known positive constant, 
and yo is a given constant. We note that if p(x) = 0, (1) becomes a Bernoulli equation. If r(x) = 0, (1) 
reduces to a first order linear ordinary differential equation. The Riccati differential equation with the case of 
p(x) = 0 or r(x) = 0 can be solved using standard methods for ordinary differential equations. In this paper, 
we consider that p(x) # 0 and r(x) # 0. Due to the important roles of Riccati differential equations in 
control systems engineering, a simple but accurate solver is desired. Providing a simple and accurate method 
for solving Riccati differential equations is the aim of this paper. 

This paper is organised as follows. We explain the problem that we want to tackle in section 2. We 
propose a multistage successive approximation method for solving Riccati differential equations in section 3. 
Results and discussion are provided in section 4. The paper is concluded in section 5. 


2. PROBLEM DESCRIPTION 

In this section, we recall an existing piecewise variational iteration method for solving Riccati 
differential equations due to Geng, ef al. [29]. The variational iteration method itself, for the general case, 
was originally proposed by He [30]-[32]. 
Considering Riccati differential (1) with initial condition (2) on domain 0 < x < X, Geng, et al. [29] took the 
correction functional 


Inr) = Yale) + Jo ACE) [PR -ORE -r ORE) — pE) az, 3) 


where Y, is a restricted variation, that is, dy, = 0; A(€) is a Lagrange multiplier, which should be determined 
optimally. Taking the variation of (3), we obtain 


Eyn) = SE) + 8 JF ACE) [EE — OnE) -r ORE — pO] dé, (4) 
Simplifying (4), we obtain 


Eyn) = Syn(x) + 5 JE ACE) [FBO] as, (5) 


Integrating (5) by parts, we have 

Syn) = 5[(1 EAA) mA] + 6 SS yn (GA (EAE, (6) 
Considering (6), we come to the following stationary conditions 

1+A(x) = 0,7'(€) = 0. (7) 
Solving stationary conditions (7), we obtain that the optimal Lagrange multiplier is 


AQ) = -1. (8) 
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Therefore, with Lagrange multiplier (8), the variational iteration method due to Geng, et al. [29] is 


In) = Ynlx) — fo [PE — OnE — ORO - PO] aé, ©) 


which is for solving Riccati differential (1) with initial condition (2) on domain 0 < x < X. 

Realising that the variational iteration method (9) produces analytical approximate solutions which 
are accurate only for points close enough to the initial position, Geng, et al. [29] implemented the method (9) 
piecewisely. By piecewisely, Geng, et al. [29] meant that the original interval I = [0, X] was subdivided into 
a finite number of subintervals I = [xj-1,)], where j = 1,2,3,...,J for a positive integer J. The width of 
each subinterval is assumed to be the same, that is, Ax = x; — x;_, for all j. With this setting, we have J + 1 
discrete points of the original interval 7 = [0, X], that is, xọ = 0, x, = Ax, x, = 2Ax, 1 Xj = JAx = X. We 
denote yjn(x) the analytical approximate solution on the jth subinterval at the nth variational iteration. 
Suppose that the maximum number of variational iterations is N, where N is a specified positive integer. This 
means that n = 0,1, 2,...,N. 

The piecewise variational iteration method due to Geng, et al. [29] works as follows. As the first 
step, on subinterval 1, = [xo, xı], we take the initialization 








Y1,0(x) = y(X%o) = y(0) = yo (10) 


and iterations 


Yas) = ine) — LP PE — ain) — Oven — PE] az, (11) 


for n = 0,1,2,...,N — 1. As shown in (10) and (11), as the second step, on subinterval I, = [x,, x2], we 
take the initialization 


Y2,0(x) = yın (%) (12) 


and iterations 


Yon) = Yan) — fe [PEE — Oyn) — rT OyFnG) - PO] aé, (13) 


for n = 0,1,2,...,N — 1. As shown in (12) and (13), as the next steps, on subintervals I; = [jaa x], where 
j =3,4,5,...,J, we take the initialization 


Yj o) = Yj-1,n Xj-1) (14) 


and iterations 


Viar) = Yin) — SL PES- a yin) -r OaE) — PO] as, (15) 


for n = 0, 1,2, ..., N — 1. Using (10)-(15), we obtain the solution on the whole domain. 

Now, the most expensive computation with the formulation of Geng, et al. [29] above lies on the 
part where we need to calculate the derivative dy, ,(¢)/dé and integrate the results from x;_, to x € J}. 
These tasks are redundant. To obtain a simpler method for solving Riccati differential equations, we need to 
modify this piecewise variational iteration method of Geng, et al. [29]. This is the problem that we aim to 
solve in this paper. 


3. PROPOSED MULTISTAGE SUCCESSIVE APPROXIMATION METHOD 

In this section, we propose a modification of the piecewise variational iteration method of Geng, et 
al. [29], so that the modified method is simpler in terms of computing implementation, yet its accuracy does 
not change. We reconsider the iterative formula of the variational iteration method (9) and do the integration 
for the term having the derivative of the unknown function, so as shown in (9) becomes 





A multistage successive approximation method for Riccati differential equations (Petrus Setyo Prabowo) 


1592 O ISSN: 2302-9285 


Yn) = yo + SORNE + rOyZ(E) + v1 dé, (16) 


where for n = 0,1,2,...,N — 1. Iterative formula (16) is Picard’s successive approximation method for 
Riccati differential (1). The iterative in (16) of the successive approximation method is equivalent to the 
iterative formula of the variational iteration method (9), but in (16) is simpler. 

Now, we shall implement the iterative (16) of the successive approximation method (SAM) 
piecewisely. The resulting method is called multistage successive approximation method (MSAM) for 
solving Riccati differential equations. Our MSAM works as follows. As the first step, on subinterval I, = 
[xo; x1], we take the initialization 





Yio) = y(xo) = y(0) = Yo (17) 


and iterations 


Viner) = Yo + N layin) + rY?) + e@] dé, (18) 


forn = 0,1,2,...,N — 1. As shown in (17) and (18), as the second step, on subinterval J, = [x1, x2], we take 
the initialization 


Y2,0(%) = yın) (19) 


and iterations 
Vons) = Yoox) + fe [a @ yen) + ryn) + pO] as, (20) 


for n = 0,1,2,...,N — 1. As shown in (19) and (20), as the second step, on subinterval I = [x-1], we 
take the initialization 


Yj o) = Yj-1,n %Gj-1) (21) 


and iterations 
Vines) = Yo) + fe, [IOyjin@ Hr Oya + PO] as, (22) 


for n = 0, 1,2, ..., N — 1. Using (17)-(22), we obtain the solution on the whole domain. 
To optimise our MSAM further, in the computer implementation, we compute the successive 
approximation formula symbolically only once: 


Yin) = Haj) +S [a Oyjn@) +r @Oyin@ + pO] as, (23) 


for n = 0, 1,2,...,N — 1. Then, the obtained symbolic formula from (23) is used to solve the Riccati 
differential equation on each subinterval J; consecutively for j = 1, 2,3, ..., J. Here, we speficy that 


Yon(Xo) = Yo (24) 


for all n. In (24) means that the given initial value is used as the starting point of solution. 


4. RESULTS AND DISCUSSION 

In this section, we provide research results on the performance tests of our proposed method and 
discuss about them. We take two computational tests, namely, a Riccati differential equation with constant 
coefficients and a Riccati differential equation involving a variable coefficient. Error on the considered 
domain is defined as the average of relative errors at all discrete points on the domain. 


4.1. Riccati differential equation with constant coefficients 
As the first test, we consider the Riccati differential equation with constant coefficients [29]: 
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OM =1+2y(x)-y?X),0<x < 4, N 


dax 
having initial condition 
y(0) =0 (26) 


The exact solution to this problem is 





1 -14+v2 
y(x) = 1 + V2tanh (v2x += log ( ae )). (27) 

Our computational experiments show that the standard SAM is not able to solve the problem on 
thewhole domain. In contrast, our proposed MSAM is able to solve the problem on the whole domain 
accurately. These phenomena are shown in Figure 1(a), where SAM and MSAM use 3 successive iterations, 
and in addition, for MSAM we use Ax = 0.1. 

To investigate further about the performance of MSAM, we record the errors and their orders of 
convergence in Tables 1-4. We obtain that the number of successive iterations in the MSAM evolution 
determines the order of convergence of the solution. One successive iteration in the MSAM evolution leads 
that MSAM is a first order method. This is observed as Ax approaches to zero, the order of convergence 
tends to one, as recorded in Table 1. Two successive iterations in the MSAM evolution leads that MSAM is a 
second order method, because as Ax approaches to zero, the order of convergence tends to two, as recorded 
in Table 2. Three successive iterations in the MSAM evolution leads that MSAM is a third order method, 
because as Ax approaches to zero, the order of convergence tends to three, as recorded in Table 3. Similarly, 
four successive iterations in the MSAM evolution leads that MSAM is a fourth order method, because as Ax 
approaches to zero, the order of convergence tends to four, as recorded in Table 4. In general, smaller Ax 
results in smaller error. Furthermore, more successive iterations results in higher order accurate method. 


Table 1. Error and order of convergence of MSAM solution for Riccati differential equation with constant 
coefficients, in which we use 1 iteration in the MSAM evolution. Error is computed on interval [0,4] 








Ax Average of Relative Error Order of Convergence 
0.25 4.420E-02 = 
0.125 2.425E-02 0.87 
0.0625 1.274E-02 0.93 
0.03125 6.536E-03 0.96 
0.015625 3.311E-03 0.98 





Table 2. Error and order of convergence of MSAM solution for Riccati differential equation with constant 
coefficients, in which we use 2 iterations in the MSAM evolution. Error is computed on interval [0,4] 








Ax Average of Relative Error Order of Convergence 
0.25 7.393E-03 = 
0.125 2.022E-03 1.87 
0.0625 5.395E-04 1.91 
0.03125 1.399E-04 1.95 
0.015625 3.565E-05 1.97 





Table 3. Error and order of convergence of MSAM solution for Riccati differential equation with constant 
coefficients, in which we use 3 iterations in the MSAM evolution. Error is computed on interval [0, 4] 








Ax Average of Relative Error Order of Convergence 
0.25 6.674E-04 - 
0.125 9.271E-05 2.85 
0.0625 1.252E-05 2.89 
0.03125 1.636E-06 2.94 
0.015625 2.094E-07 2.97 





Table 4. Error and order of convergence of MSAM solution for Riccati differential equation with constant 
coefficients, in which we use 4 iterations in the MSAM evolution. Error is computed on interval [0,4] 








Ax Average of Relative Error Order of Convergence 
0.25 7.349E-05 = 
0.125 5.049E-06 3.86 
0.0625 3.388E-07 3.90 
0.03125 2.205E-08 3.94 
0.015625 1.407E-09 3.97 
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4.2. Riccati differential equation involving a variable coefficient 
As the second test, let us consider the Riccati differential equation with a variable coefficient [29]: 


dy(x) _ 





7 =1+x7-y7(x),0<x <4, (28) 
with initial condition 
y(0)=1 (29) 
The exact solution to this problem is 
ex” 
yx) =xt (30) 


1+ fo e7 $2 dg 


For this second test, the standard SAM is not able to solve the problem, even for the first half [0,2] 
of the given domain [0,4], as shown in Figure 1(b). In contrast, our proposed MSAM is able to solve the 
problem on the whole domain [0,4] accurately. For this Figure 1(b), SAM and MSAM use 3 successive 
iterations, and in addition, for MSAM we take Ax = 0.1. 

To investigate further about the performance of MSAM in solving this second test case, we record 
the errors and their orders of convergence in Tables 5-8. The behaviour of MSAM in this test case is 
consistent with that of MSAM in the previous test case. We obtain that the number of successive iterations in 
the MSAM evolution determines the order of convergence of the solution. One successive iteration in the 
MSAM evolution leads that MSAM is a first order method. This is observed as Ax approaches to zero, the 
order of convergence tends to one, as recorded in Table 5. Two successive iterations in the MSAM evolution 
leads that MSAM is a second order method. This is observed as Ax approaches to zero, the order of 
convergence tends to two, as recorded in Table 6. Three successive iterations in the MSAM evolution leads 
that MSAM is a third order method. This is observed as Ax approaches to zero, the order of convergence 
tends to three, as recorded in Table 7. Similarly, four successive iterations in the MSAM evolution leads that 
MSAM is a fourth order method, because as Ax approaches to zero, the order of convergence tends to four, 
as recorded in Table 8. Again, in general, smaller Ax results in a more accurate method. Furthermore, more 
number of successive iterations makes MSAM to be higher order accurate. 

With our accurate results in this paper, we are confident that MSAM is a reliable method to be used 
for other kinds of initial value problems. The idea of MSAM could be adapted for continuous versions of 
discrete problems in computing [33], and it could be extended to be analysed using advanced mathematical 
tools, such as, algebraic geometry and conformal mapping. 


4 4 - n 
-Exact solution j a 
2 eee iin eet gl LesaMeontion y ae 
_ we ae __ | Z MSAM solution | a 
—Exact solution a4 yet ar | 
-2” SAM solution Nes. P| qpe | 
~MSAM solution ~ 
-4 o 
0 1 2 3 4 0 1 2 3 4 
X x 
(a) (b) 


Figure 1. Exact, SAM, and MSAM solutions on interval [0, 4]. SAM solution is accurate only at points close 
to the initial condition. MSAM solution coincides graphically with the exact solution, (a) Results of the first 
test case, (b) Results of the second test case 


Table 5. Error and order of convergence of MSAM solution for Riccati differential equation involving a 
variable coefficient, in which we use 1 iteration in the MSAM evolution on interval [0,4] 








Ax Average of Relative Error Order of Convergence 
0.25 3.544E-02 - 
0.125 1.764E-02 1.01 
0.0625 8.802E-03 1.00 
0.03125 4.396E-03 1.00 
0.015625 2.197E-03 1.00 
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Table 6. Error and order of convergence of MSAM solution for Riccati differential equation involving a 
variable coefficient, in which we use 2 iterations in the MSAM evolution on interval [0,4] 








Ax Average of Relative Error Order of Convergence 
0.25 3.408E-02 - 
0.125 4.625E-03 2.88 
0.0625 9.436E-04 2.29 
0.03125 2.165E-04 2.12 
0.015625 5.202E-05 2.06 





Table 7. Error and order of convergence of MSAM solution for Riccati differential equation involving a 
variable coefficient, in which we use 3 iterations in the MSAM evolution on interval [0,4] 








Ax Average of Relative Error Order of Convergence 
0.25 6.406E-03 - 
0.125 6.894E-04 3.22 
0.0625 7.439E-05 3.21 
0.03125 8.555E-06 3.12 
0.015625 1.024E-06 3.06 





Table 8. Error and order of convergence of MSAM solution for Riccati differential equation involving a 
variable coefficient, in which we use 4 iterations in the MSAM evolution on interval [0,4] 








Ax Average of Relative Error Order of Convergence 
0.25 2.393E-03 - 
0.125 1.019E-04 4.55 
0.0625 5.294E-06 4.27 
0.03125 3.007E-07 4.14 
0.015625 1.790E-08 4.07 





5. CONCLUSION 

We have proposed successfully a new multistage successive approximation method for solving 
Riccati differential equations. The main advantage of the proposed method is that it is simpler than the 
existing variational iteration method for solving Riccati differential equations. The proposed method is 
analytically equivalent to the existing method, but simpler in terms of computing implementation. We have 
tested the performance of the multistage successive approximation method in cases that their exact solutions 
are known. The order of accuracy of the proposed method can be made higher by simply taking more number 
of successive iterations in the multistage evolution. Obviously, taking smaller size of each subinterval and 
taking more number of iterations in the multistage evolution lead that our proposed method produces smaller 
error and it becomes higher order accurate. With these results, in cases that the exact solutions to Riccati 
differential equations are not known, we are confident to propose the use of the multistage successive 
approximation method for solving them to obtain accurate approximate solutions. 
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